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Abstract 

Wc study the axisymmetric propagation of a viscous gravity current over a deep porous 
medium into which it also drains. A model for the propagation and drainage of the current 
is developed and solved numerically in the case of constant input from a point source. In this 
case, a steady state is possible in which drainage balances the input, and we present analytical 
expressions for the resulting steady profile and radial extent. We demonstrate good agreement 
between our experiments, which use a bed of vertically aligned tubes as the porous medium, 
and the theoretically predicted evolution and steady state. However, analogous experiments 
using glass beads as the porous medium exhibit a variety of unexpected behaviours, including 
overshoot of the steady-state radius and subsequent retreat, thus highlighting the importance 
of the porous medium geometry and permeability structure in these systems. 



1 Introduction 

Gravity currents are primarily horizontal fluid flows driven by a density difference between the 
intruding and ambient fluids. These flows are common in natural systems and industrial processes 
and describe, for example, the spread of cold air into a room, the dispersal o f pollutants from 
an industrial spill, and the flow of snow and debris avalanches ( Huppert . 20061 ). Many previous 
studies have examined in detail the propagation of currents along impermeable boundaries; here 
we consider flow over porous substrates through which these currents can also drain. 

Two-dimensional gravity currents propagating over porous media have been addressed both 
theoretically and experimentally by several authors. For curre nts flowing over thin porous sub- 
strates, only the wei g ht of the overlying fluid drives drainage (Thomas. Marino fc Linden . 19981 : 
Ungarish fc Huppertl . l200d : iMarino fc Thomasl . 120021 : IPritchard. Woods Hoed. l200lh. In con- 



trast, for gravity currents propagating over deep porous media, lActon. Huppert fe Worsterl (|200ll ) 
showed that both the hydrostatic pressure of the fluid in the current and the weight of the fluid 
within the porous medium drive drainage. They used this description of drainage in a model of 
experiments in which low Reynolds number gr avity currents spread over a deep porous layer in two 
dimensions. iThomas. Marino &: Linden! (|2004l ) used this drainage law to describ e their experiments 
on the propagation of high Reynolds number currents over deep porous media. IPritchard &; Hogg 
( 20021 ) have also applied the same drainage law to their examination of gravity currents propagating 
within a porous medium overlying a deep layer of lower permeability. Simil ar studies have also been 
cond ucted that consider two-phase flow within the porous medium (e.g. iHussein. Jin &: Weaverl . 
2002J), but these effects are beyond the scope of the present study. 
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Axisymmetric gravity currents propagating over porous media have been studied primarily as 
microscale flows in which capill ary forces drive the drainage and therefore the wetting prop erties of 
the medium are important ( e.g. 



ry iorces drive tne dra inage and tncr ciore tne wetting prop erties 01 
Davis fc Hockind . ll99Sl . 120001 : iKumar fc DeshDandd . 120061 ) . At the 



macroscale, iPritchard et al\ (|200ll ) considered gravity-driven drainage of an axisymmetric current 
flowing through a porous medium overlying a thin layer of lower permeability. In both geometries, 
previous experiments only involved currents of fixed volume, whereas our experiments explore the 
fixed flux case. 

Here we examine the axisymmetric propagation of a macroscopic viscous gravity current over a 
deep porous medium. Ou r model uses lubrication theory for flow within the current, the drainage 
law of I Acton et~aH (|200lh . and Darcy flow within the porous medium. While the full spatial and 



temporal evolution of the current can only be obtained numerically, an analytical expression for the 
steady-state extent and profile of a current fed by a constant input of fluid is found. Additionally, 
we develop scaling laws describing the propagation of the current. Our experimental setup in 
which a gravity current fed by a constant flux of golden syrup spreads across a bed of vertically 
aligned straws conforms closely to the assumptions of our model, so the scaling laws provide a 
good collapse of all data onto a curve in agreement with the numerical solution. In contrast to 
these well-behaved currents, we describe the non-ideal behaviour observe d in experiments using 
glycerin and glass beads similar to the system used by lActon et al\ (120011 ) . We propose that the 
axisymmetric geometry makes the currents particularly sensitive to any non-uniformities of the 
porous medium which leads to the disagreement between these experiments and the theory. 



2 Theoretical model 

We consider the axisymmetric spreading of a fluid of kinematic viscosity v and density p into an 
ambient fluid of density p a « p and viscosity v a « v. As shown schematically in figured! fluid is 
supplied at the origin and spreads radially over a porous medium with porosity <j) and permeability 
k into which the fluid drains. We consider the general case in which the volume of fluid increases as 
qt a , where t is time and q and a are constants. After a brief initial stage, the radial extent of the 
current rw(t) is much greater than its height h(r,t), and in this limit we apply the approximations 
of lubrication theory: velocity within the current is assumed to be predominantly horizontal and 
pressure within the current is assumed to be hydrostatic. Under this approximation, viscous flow 
within the current is driven by radial gradients of its thickness. We apply conditions of no slip at 
z = and no tangential stress at z = h to determine the horizontal fluid flux q^ = — (g/3u) rh s dh/dr. 
Conservation of fluid mass through an infinitesimal control volume of the gravity current gives the 
equation 

dh q 1 d I , -idh\ , n N 
at Sv r or \ or ) 

governing the current's structure and evolution, where w(r, 0,t) is the drainag e veloc ity from the 



base of the current into the underlying porous medium. Following lActon et al\ ()200ll ). we assume 
that drainage into the porous medium is driven both by the weight of the draining fluid and the 
hydrostatic pressure of the fluid within the current, giving 

qk I h\ dl 
w ( rA t) = + = -</>- (2) 

where l(r,t) is the depth of the fluid within the porous medium. 

In this analysis, we have made a few assumptions that merit further examination. First, our 
assumption of no slip at the porous medium surface is valid when the current height is much greater 
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Figure 1: An illustration of the theoretical and experimental geometry considered for the axisym- 
metric spreading of a viscous gravity current over a porous medium. 



than the pore size because the presence of a slip veloc ity is equivalent to extending the fluid region 
a dis tance of less than one pore size into the medium (jBeavers fe Joseph! . 119671 ; iLe Bars Worsted . 
20061 ). Near the nose of the current or for currents flowing over very rough substrates apparent slip 
may be importa nt. Secondly , we ha ve assumed that surface tension is negligible in the drainage law. 
As discussed in lActon et al\ (|200ll ) , this is accurate as long as the pressure due to surface tension 
is much less than the hydrostatic pressure; equivalently the capillary rise height h c ~ ^y/pga, where 
7 is the surface tension and a is the pore radius, must be much larger than h. For our experimental 
setup h c « 2 mm, much less than typical current heights. Finally, equation ((2J) assumes that flow 
within the porous medium is single-phase and that the porous matrix is stationary; therefore we 
can ignore flow of the displaced fluid and assume that k is constant in time. This is an accurate 
assumption when the wetting properties of the displaced and displacing fluids with respect to the 
porous matrix are similar and when the displaced fluid is inviscid. Equation ([2]) also assumes 
that fluid flow is predominantly vertical, implicitly neglecting the potential for the Rayleigh-Taylor 
instability. 

The governing equations (JT]) and ([2]) are subject to one boundary condition specifying the flux 
near the origin and another requiring zero flux through the nose of the current rj\r(i). Respectively, 
they are 



lim 

r-»0 



2irr 



gh dh 
3v dr 



= -aqt 



a-l 



and 



2-Kr 



gh 3 dh 
2>u dr 



0. 



(3a, b) 



Jr N 



We note that boundary condition (|3]a) along with the evolution equation ([I]) is equivalent to a 
statement of global mass conservation, namely 



rr N rt rr N 

qt a = 2TT / rhdr-2-K / rw(r,0,t) dr dt 
Jo Jo Jo 



(4) 



We non-dimensionalize equations (HI)-© by introducing horizontal, vertical and temporal scales 
Sh, Sy and St given by 



(a-l)/2(a-5) 



V \3g^ 



l/(a-5) 



and St 



3qv) 



l/(«-5) 



(5a, b, c) 



where V = gk\v is the characteristic drainage velocity in the porous medium. By introducing 
dimensionless variables 



H-h/Sy, L-l/Sy, R-rjSn and T-t/Sx, 



(6a, b, c, d) 
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the equations governing the dimensionless height H(R,T) and depth L(R,T) of the intruding fluid 
become 



dH 1 d 
~df~RdR 



(**tf)~K) 



and 



V dT \ L J 
These are subject to the scaled boundary conditions 



lim 



dR 



-aT a ~ l and 



2irRH d 



~dR 



(7) 
(8) 

(9a, b) 



For a fixed flux at the origin (a = 1) the scalings simplify to 



1/2 



i/i 



and /St = f -=-, . I 



5\l/4 



(10a, b, c) 



This system admits a steady state in which drainage exactly balances the material input. In the 
long-time limit, the depth of drainage greatly exceeds the height of the current L » H and the 
current has the steady profile 



H=[R 2 



R% 



2R 2 N In (R/R N )] 



i/i 



(11) 



plott ed in figure [2j The logarithmic singularity at R = accounts for the finite flux there (jHuppertl . 



1982 ), Balancing the external flux with the drainage flux, we find the steady-state extent 



.R/v(t -*■ oo) = 7T 



-1/2 



0.564. 



(12) 



We note that Pritchard et al. determined an analytical solution for the steady state of 

currents in a similar system. They considered a two-dimensional current fed by a constant source 
at the origin spreading through a porous medium of high permeability and underlain by a very 
thin low permeability porous layer. Comparing their figure 3 and equation (2.20) with our figure 
[5] and equation (|11|) highlights the importance of the specific system to the current behaviour. 
Whereas in our system, the steady-state profile of the current has an inflection point where the 
curvature switches from being positive near the source to negative near the nose, in their system 
the steady-state current surface is concave upwards near the nose. 



3 Numerical solution 

The full time evolution of the current height H (R,T) and depth L(R,T) is found numerically 
by integrating equations ([7]) and (jHJ) on a uniform grid with spacing 0.001. We compute the new 
drainage depth L n+ \ = L(R,T n+ i) from equation (jHJ) using the height H n = H (R,T n ) and depth 
L n = L(R,T n ) from the previous time step. We then use L n+ \ to compute the drainage velocity 
on the right-hand side of equation (J7]) and solve for the height from equation ([7]) using the control 
volu me (or flux conse rvative) method in space and a Crank-Nicholson (semi- imp licit) scheme in 
time ( Patankarl . 1980l ). In this computation, (H n ) 3 is our initial estimate for if 3 in the non-linear 



term on the left-hand side of equation ([?]). We update the height as described above and use the 
new estimate for H in the non-linear term, iterating this process until the updated value of H 
converges. This converged value is H n+ \. Finally, we proceed to the next time step. 
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Figure 2: The steady-state profile of dimensionless height H versus dimensionless radius R. 



We introduce fluid into the current by assigning a constant flux at the left-hand boundary of 
the first control volume {R = 0) in the discretized equations. The right-hand boundary of our grid 
is impermeable to fluid flow and is positioned beyond the steady-state extent. The initial condition 
is an empty box with no fluid. We record the nose of the current as the position where H(R) 
falls below a prescribed small tolerance: the height beyond this point is set to 0. This condition 
is necessary because the drainage velocity is ill-defined when L = and H ^ 0, as occurs near the 
current nose at the beginning of a time step. 

We have tested the sensitivity of the numerical solution to the choice of grid spacing, time-step 
size and height tolerance, and found the results to be relatively insensitive to these parame ters. Ad- 
dition ally, we have tested our numerical results for non-draining currents against those of iHuppert 
(119821 ) . and find good agreement for both the constant volume and constant flux cases. 

Figure [3] shows three calculated profiles (curves) at different times for a numerically simulated 
current with porosity <fi = 0.907, similar to that of our experimental setup. We note that, as the 
system approaches steady state, the volume of fluid residing within the porous medium becomes 
much larger than the volume in the current above the medium. The solid curve in figure [5] shows 
how the extent increases with time, approaching the steady-state value predicted by equation (fl2]) , 



4 Experiments 

We performed a series of experiments using Lyle's golden syrup as the viscous intruding fluid and 
a bed of vertically oriented drinking straws of radius r s = 0.29 ± 0.01 cm as the underlying porous 
medium, as shown in figure [H Lyle's golden syrup was used as the working fluid because its 
viscosity of v > 400 cm 2 s _1 (as measured by a U-tube viscometer) results in currents with heights 
much greater than the surface topography of the porous medium. The simple geometry of the 
porous medium ensures strictly vertical drainage flow and allows for a comparison between the 
experimentally measured and the theoretically predicted permeability. 

We measured the permeability of the porous medium by conducting drainage experiments in 
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Figure 3: Dimensionless height H and depth L profiles from the numerical solution with <f> = 0.907 
(curves) and experiment 9 (symbols) at dimensionless times T = 0.216 (solid curve and □), 0.592 
(dashed curve and o), and 1.72 (dotted curve and 0)- Data for R < 0.1 are not plotted as effects 
from the finite width and coiling instability of the fluid source are most pronounced there. See 
discussion of experimental errors in §4 




Figure 4: Lyle's golden syrup spreading out over and draining into a bed of drinking straws. 
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3riment 


q (cm 3 s 1 ) 


v (cm 2 s 1 ) 


Stj (cm) 


Sy (cm) 


St (s) 


1 


1.06 ±0.01 


453 


8.78 


1.10 


80.0 


2 


4.18 ±0.03 


453 


17.4 


1.55 


113 


3 


2.09 ±0.01 


453 


12.3 


1.30 


94.7 


4 


9.83 ±0.01 


453 


26.7 


1.92 


139 


5 


6.74 ±0.01 


453 


22.1 


1.75 


127 


6 


2.31 ±0.01 


401 


12.2 


1.30 


83.4 


7 


1 27 ± 01 


401 


9.05 


1.12 


71.9 


8 


2.00 ±0.01 


401 


11.3 


1.25 


80.4 


9 


6.11 ±0.01 


401 


19.8 


1.65 


106 


10 


7.37 ±0.01 


401 


21.8 


1.73 


111 


11 


6.63 ±0.01 


401 


20.7 


1.69 


109 


12 


20.35 ± 0.03 


401 


36.2 


2.24 


144 



Table 1: Summary of the experimental parameters. For each experiment, the permeability was 
assumed to be k - 6.36 ± 0.04 x 10~ 3 cm 2 . The uncertainty in viscosity is described in the text. 

which syrup with v = 453 cm 2 s _1 and p = 1.5 g cm~ 3 was maintained at a constant height h = 10 cm 
above the porous medium within a large cylinder of radius rc = 5.75 cm as it drained through the 
straws of length / = 20 cm. By measuring the mass flux dM/dt through the straws with a digital 
scale connected to a computer, we obtained the drainage velocity 

v \ I I P^Tq 

and consequently a measure of the permeability of the porous medium 

k exp = 6.36 ± 0.04 x 10~ 3 cm 2 . (14) 

The uncertainty in this value comes from estimating dM/dt from the measured mass versus time, 
which has an uncertainty of ±0.02 g s _1 . This experime ntal v a lue c an be compared with the 



theoretical permeability for aligned capillary tubes given by iBearl (|1972l ) 



6r 2 

/c = = 9.5 ± 0.3 x 10~ 3 cm 2 , (15) 

where <f> = ir\/3/6 - 0.907 is the packing fraction of the straws for hexagonal close packing. We 
attribute the approximately 30% discrepancy between the measured and theoretical values to a slow 
leakage of golden syrup through the interstices between the straws and imperfections in the straw 
packing that produced a porosity not equal to that of a close packing. In the following analysis of 
the experimental data, we use the measured permeability. 

For each experiment, a fixed flux of syrup was supplied at the origin from a reservoir maintained 
at a constant gravitational head. The mass flux was measured with a digital balance connected 
to a computer prior to the initiation of each experiment. The flux, viscosity and resultant scaling 
laws are summarised in tabled] for each experiment. 

Digital images of the side profile of each experiment were made at regular intervals (see figure S]), 
and later analysed to obtain the radial extent and height profiles of each current. A comparison 
between the scaled radial extent of each current and the numerical solution to equations (JT])-® is 
shown in figure The dotted curves represent a ±10% error bound in Sh applied to the numerical 



7 




Figure 5: The dimensionless radial extent Rn versus dimensionless time T for the experiments 
listed in tableland our numerical solution (solid curve). The dotted curves show the ±10% error 
bounds. The symbols correspond to experiments: 1 (□), 2 (o), 3 (A), 4 (o), 5 (v), 6 (+), 7 (x), 8 
(■), 9 (.), 10 (A), 11 (♦), and 12 (▼). 

extent (solid curve). The error in St is not represented in the plot. Uncertainty in v and k are the 
main contributors to the overall uncertainty, as the error in q is less than 1%. Although the viscosity 
was measured regularly throughout the set of experiments, the large range of values obtained and 
the known large temperature dependence of the viscosity (20% per °C) result in a large uncertainty 
in the actual viscosity. Due to this uncertainty in the viscosity and the discrepancy between the 
theoretical and measured permeabilities, we estimate the total error to be about ±10%. Within 
the error bounds, the collapse of the scaled data and agreement with the numerical solution is 
good. We also obtained height profiles from the images of experiment 9, which are compared to 
the numerical profiles at the same scaled times in figure [3j For clarity error bounds are not plotted, 
but the uncertainty is again ±10%. Although the finite width and coiling instability of the fluid 
source cause some discrepancy, the overall agreement is good. 



5 Discussion 

We have shown that for a fixed flux of golden syrup flowing across a bed of vertically ali gned straws a 



simple model based upon lubrication theory and the drainage law of lActon et al\ (|200ll ) can describe 
the current propagation and steady state. In contrast, experiments conducted using glycerin as the 
working fluid and ~ 3 mm diameter spherical glass beads as the porous medium (detailed results not 
included here) exhibited non-ideal behaviour that violated a number of assumptions in our model. 
In particular, most currents had a scalloped front as the current propagated across the beads (figure 
Eh), complicating measurement of the current radius. Many were also non-axisymmetric as shown 
in figure [6]&. Finally, all of the glycerin currents exhibited a maximum extent from which the current 
nose then retreated (figure [6]c) . 

We attribute these non-ideal behaviours primarily to a sensitive dependence on the character- 
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Figure 6: Examples of non-ideal behaviour from three gravity currents using glycerin (dark) input 
at the upper right corner and spreading across beads (light); (a) large and small scallops at the 
leading edge; (b) non-axisymmetric spreading (the quarter circle is equiaxial); and (c) a region of 
beads over which the current advanced and subsequently receded. 

istics and geometry of the underlying porous medium. For example, at the nose of the current the 
thickness is small and may be comparable to the surface topography of the porous medium. This 
could cause the front to stick on surface asperities, producing the scalloped edge. Additionally, 
the currents are sensitive to inhomogeneities in the bead packing, and thus the permeability of the 
medium, due to the strong influence of permeability on the drainage velocity in equation ([2]) and 
on the scaling laws in equation (1101) . This may have contributed to the scalloped front and the 
non-axisymmetric propagation. The non-axisymmetry also could have arisen from a bead surface 
that was not sufficiently level. Although care was taken to level the surface, we cannot exclude this 
possibility. These hypotheses are supported by our experiments using golden syrup and straws, a 
level and uniform porous medium, and could be tested by conducting more experiments using, for 
example, glycerin and straws or smaller beads. 

For the roll-back phenomenon, we have no simple explanation. However, we can rule out some 
possibilities. First, we verified that horizontal flow within and immediately above the surface of the 
porous medium was negligible as assumed in our model (no-slip condition). Powdered dye placed in 
several small piles on the bead surface along the path of the current was picked up by the draining 
fluid and carried purely vertically into the beads. Secondly, the geometry of the glycerin experiments 
afforded us a cross-sectional view of the current and draining fluid from which we observed a uniform 
saturation of the beads. This supports our assumption of a constant permeability, though as we 
could not observe the interior of the porous medium, we cannot completely rule out these effects in 
the bulk of the flow. Because the roll-back phenomenon was not observed in the golden syrup and 
straws system, we think that it is related to the specific combination of fluid and porous medium 
properties. Again, this could be tested with experiments involving different fluids and porous media. 
Finally, we note that no experimental evidence for a Rayleigh-Taylor instability was found at the 
lower interface of the current on the time scales over which the experiments were conducted. This 
observation implies that, at least here, vertical drainage is the dominant factor controlling radial 
spreading of the current. 



The contrast between our experiments using glycerin and beads and those of lActon et all (|200ll ) 
using glycerin and beads in a linear geometry with a fixed fluid volume suggests that some char- 
acteristic of either the axisymmetric geometry or the fixed fluid input results in currents that are 
much more sensitive to the properties of the porous medium. For example, axisymmetric currents 
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have a much longer front and therefore a larger nose area than linear currents. Therefore their 
spreading is more strongly influenced by surface roughness and the failure of our model assump- 
tions near the nose. To explore these ideas further, we suggest conducting fixed flux experiments 
in the linear geometry and fixed volume experiments in the axisymmetric geometry using different 
porous media. 

The sensitive dependence of propagation and drainage of the current on the spatial structure of 
the permeability and the surface topography suggests that further studies are needed to characterise 
fluid flow in these situations. Nonetheless, our model provides a simple framework to estimate the 
evolution of the current over time and the maximum extent at steady state for currents flowing 
over simple porous media. 
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periments. This research was partially supported by the U.S. Department of Energy (DE-FG02- 
05ER15741). MJS acknowledges support from a U.S. National Science Foundation Graduate Re- 
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Royal Institute of Technology, and NORDITA in Stockholm. 
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